﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using SmartMathLibrary.Blas.Engine;

namespace SmartMathLibrary.Blas
{
    /// <summary>
    /// This class provides an interface to the Blas level 3 routines.
    /// </summary>
    [Serializable]
    public static class ManagedBlasLevel3
    {
        /// <summary>
        /// Sum of the product of two general matrices and a scaled matrix. 
        /// </summary>
        /// <param name="transA">Specifies the operation to be performed on the matrix a.</param>
        /// <param name="transB">Specifies the operation to be performed on the matrix b.</param>
        /// <param name="m">The number of rows in the matrix a transformed as specified by transA, and the matrix c.</param>
        /// <param name="n">The number of columns in the matrix b transformed as specified by transB, and the matrix c.</param>
        /// <param name="k">The number of columns in the matrix a transformed as specified by transA, and the number of rows of the matrix b transformed as specified by transB.</param>
        /// <param name="alpha">The scalar used to multiply the matrix-vector product.</param>
        /// <param name="a">Reference to the first element in a one-dimensional array containing the elements of the first matrix.</param>
        /// <param name="lda">Leading dimension of a.</param>
        /// <param name="b">Reference to the first element in a one-dimensional array containing the elements of the second matrix.</param>
        /// <param name="ldb">Leading dimension of b.</param>
        /// <param name="beta">The scalar used to multiply y.</param>
        /// <param name="c">Reference to the first element in a one-dimensional array containing the elements of the third matrix.</param>
        /// <param name="ldc">Leading dimension of c.</param>
        public static unsafe void Dgemm(TransposeOperation transA, TransposeOperation transB, int m, int n, int k,
                                        double alpha,
                                        double[] a, int lda, double[] b, int ldb, double beta, double[] c, int ldc)
        {
            fixed (double* pA = a)
            {
                fixed (double* pB = b)
                {
                    fixed (double* pC = c)
                    {
                        sbyte* pTransA;

                        if ((transA == TransposeOperation.None) || (transA == TransposeOperation.Conjugate))
                        {
                            sbyte tempuri = (sbyte) 'N';
                            pTransA = &tempuri;
                        }
                        else
                        {
                            sbyte tempuri = (sbyte) 'T';
                            pTransA = &tempuri;
                        }

                        sbyte* pTransB;

                        if ((transB == TransposeOperation.None) || (transB == TransposeOperation.Conjugate))
                        {
                            sbyte tempuri = (sbyte) 'N';
                            pTransB = &tempuri;
                        }
                        else
                        {
                            sbyte tempuri = (sbyte) 'T';
                            pTransB = &tempuri;
                        }

                        ManagedBlasProvider.dgemm_(pTransA, pTransB, &m, &n, &k, &alpha, pA, &lda, pB, &ldb, &beta, pC,
                                                   &ldc);
                    }
                }
            }
        }

        /// <summary>
        /// Sum of the product of a symmetric and a general matrix and a scaled matrix. 
        /// </summary>
        /// <param name="side">Specifies on which side the symmetric matrix a is to be multiplied.</param>
        /// <param name="triangleMode">Specifies whether the elements of the matrix a are stored in the 
        /// upper or lower triangular part.</param>
        /// <param name="m">The number of rows in the matrix a transformed as specified by transA, and the matrix c.</param>
        /// <param name="n">The number of columns in the matrix b transformed as specified by transB, and the matrix c.</param>
        /// <param name="alpha">The scalar used to multiply the matrix-vector product.</param>
        /// <param name="a">Reference to the first element in a one-dimensional array containing the elements of 
        /// the first matrix.</param>
        /// <param name="lda">Leading dimension of a.</param>
        /// <param name="b">Reference to the first element in a one-dimensional array containing the elements 
        /// of the second matrix.</param>
        /// <param name="ldb">Leading dimension of b.</param>
        /// <param name="beta">The scalar used to multiply y.</param>
        /// <param name="c">Reference to the first element in a one-dimensional array containing the elements 
        /// of the third matrix.</param>
        /// <param name="ldc">Leading dimension of c.</param>
        public static unsafe void Dsymm(MatrixOperationSide side, MatrixTriangleMode triangleMode, int m, int n,
                                        double alpha,
                                        double[] a, int lda, double[] b, int ldb, double beta, double[] c, int ldc)
        {
            fixed (double* pA = a)
            {
                fixed (double* pB = b)
                {
                    fixed (double* pC = c)
                    {
                        sbyte* pSide;

                        if (side == MatrixOperationSide.Left)
                        {
                            sbyte tempuri = (sbyte) 'L';
                            pSide = &tempuri;
                        }
                        else
                        {
                            sbyte tempuri = (sbyte) 'R';
                            pSide = &tempuri;
                        }

                        sbyte* mode;

                        if (triangleMode == MatrixTriangleMode.Lower)
                        {
                            sbyte tempuri = (sbyte) 'L';
                            mode = &tempuri;
                        }
                        else
                        {
                            sbyte tempuri = (sbyte) 'U';
                            mode = &tempuri;
                        }

                        ManagedBlasProvider.dsymm_(pSide, mode, &m, &n, &alpha, pA, &lda, pB, &ldb, &beta, pC, &ldc);
                    }
                }
            }
        }

        /// <summary>
        /// Performs a rank 2k update of a symmetric matrix. 
        /// </summary>
        /// <param name="triangleMode">Specifies whether the elements of the matrix a are stored in the upper or 
        /// lower triangular part.</param>
        /// <param name="trans">Specifies the operation to be performed on the matrix a.</param>
        /// <param name="n">The number of rows and columns in the matrix c.</param>
        /// <param name="k">The number of columns in the matrix a transformed as specified by transA.</param>
        /// <param name="alpha">The scalar used to multiply the matrix-matrix product.</param>
        /// <param name="a">Reference to the first element in a one-dimensional array containing the elements 
        /// of the first matrix.</param>
        /// <param name="lda">Leading dimension of a.</param>
        /// <param name="b">Reference to the first element in a one-dimensional array containing the elements of 
        /// the second matrix.</param>
        /// <param name="ldb">Leading dimension of b.</param>
        /// <param name="beta">The scalar used to multiply y.</param>
        /// <param name="c">Reference to the first element in a one-dimensional array containing the 
        /// elements of the third matrix.</param>
        /// <param name="ldc">Leading dimension of c.</param>
        public static unsafe void Dsyr2k(MatrixTriangleMode triangleMode, TransposeOperation trans, int n, int k,
                                         double alpha,
                                         double[] a, int lda, double[] b, int ldb, double beta, double[] c, int ldc)
        {
            fixed (double* pA = a)
            {
                fixed (double* pB = b)
                {
                    fixed (double* pC = c)
                    {
                        sbyte* pTrans;

                        if ((trans == TransposeOperation.None) || (trans == TransposeOperation.Conjugate))
                        {
                            sbyte tempuri = (sbyte) 'N';
                            pTrans = &tempuri;
                        }
                        else
                        {
                            sbyte tempuri = (sbyte) 'T';
                            pTrans = &tempuri;
                        }

                        sbyte* mode;

                        if (triangleMode == MatrixTriangleMode.Lower)
                        {
                            sbyte tempuri = (sbyte) 'L';
                            mode = &tempuri;
                        }
                        else
                        {
                            sbyte tempuri = (sbyte) 'U';
                            mode = &tempuri;
                        }

                        ManagedBlasProvider.dsyr2k_(mode, pTrans, &n, &k, &alpha, pA, &lda, pB, &ldb, &beta, pC, &ldc);
                    }
                }
            }
        }

        /// <summary>
        /// Performs a rank k update of a symmetric matrix. 
        /// </summary>
        /// <param name="triangleMode">Specifies whether the elements of the matrix a are stored in the upper or 
        /// lower triangular part.</param>
        /// <param name="trans">Specifies the operation to be performed on the matrix a.</param>
        /// <param name="n">The number of rows and columns in the matrix c.</param>
        /// <param name="k">The number of columns in the matrix a transformed as specified by transA.</param>
        /// <param name="alpha">The scalar used to multiply the matrix-matrix product.</param>
        /// <param name="a">Reference to the first element in a one-dimensional array containing the elements 
        /// of the first matrix.</param>
        /// <param name="lda">Leading dimension of a.</param>
        /// <param name="beta">The scalar used to multiply y.</param>
        /// <param name="c">Reference to the first element in a one-dimensional array containing the elements 
        /// of the third matrix.</param>
        /// <param name="ldc">Leading dimension of c.</param>
        public static unsafe void Dsyrk(MatrixTriangleMode triangleMode, TransposeOperation trans, int n, int k,
                                        double alpha,
                                        double[] a, int lda, double beta, double[] c, int ldc)
        {
            fixed (double* pA = a)
            {
                fixed (double* pC = c)
                {
                    sbyte* pTrans;

                    if ((trans == TransposeOperation.None) || (trans == TransposeOperation.Conjugate))
                    {
                        sbyte tempuri = (sbyte) 'N';
                        pTrans = &tempuri;
                    }
                    else
                    {
                        sbyte tempuri = (sbyte) 'T';
                        pTrans = &tempuri;
                    }

                    sbyte* mode;

                    if (triangleMode == MatrixTriangleMode.Lower)
                    {
                        sbyte tempuri = (sbyte) 'L';
                        mode = &tempuri;
                    }
                    else
                    {
                        sbyte tempuri = (sbyte) 'U';
                        mode = &tempuri;
                    }

                    ManagedBlasProvider.dsyrk_(mode, pTrans, &n, &k, &alpha, pA, &lda, &beta, pC, &ldc);
                }
            }
        }

        /// <summary>
        /// Product of a triangular and a general matrix. 
        /// </summary>
        /// <param name="side">Specifies on which side the triangular matrix a is to be multiplied.</param>
        /// <param name="triangleMode">Specifies whether the elements of the matrix a are stored in the upper or lower triangular part.</param>
        /// <param name="transA">Specifies the operation to be performed on the matrix a.</param>
        /// <param name="diag">Specifies whether or not a is unit triangular.</param>
        /// <param name="m">The number of rows in the matrix a transformed as specified by transA, and the matrix c.</param>
        /// <param name="n">The number of columns in the matrix b transformed as specified by transB, and the matrix c.</param>
        /// <param name="alpha">The scalar used to multiply the matrix-vector product.</param>
        /// <param name="a">Reference to the first element in a one-dimensional array containing the elements of the first matrix.</param>
        /// <param name="lda">Leading dimension of a.</param>
        /// <param name="b">Reference to the first element in a one-dimensional array containing the elements of the second matrix.</param>
        /// <param name="ldb">Leading dimension of b.</param>
        public static unsafe void Dtrmm(MatrixOperationSide side, MatrixTriangleMode triangleMode,
                                        TransposeOperation transA,
                                        MatrixDiagonalMode diag, int m, int n, double alpha, double[] a, int lda,
                                        double[] b, int ldb)
        {
            fixed (double* pA = a)
            {
                fixed (double* pB = b)
                {
                    sbyte* pTrans;

                    if ((transA == TransposeOperation.None) || (transA == TransposeOperation.Conjugate))
                    {
                        sbyte tempuri = (sbyte) 'N';
                        pTrans = &tempuri;
                    }
                    else
                    {
                        sbyte tempuri = (sbyte) 'T';
                        pTrans = &tempuri;
                    }

                    sbyte* mode;

                    if (triangleMode == MatrixTriangleMode.Lower)
                    {
                        sbyte tempuri = (sbyte) 'L';
                        mode = &tempuri;
                    }
                    else
                    {
                        sbyte tempuri = (sbyte) 'U';
                        mode = &tempuri;
                    }

                    sbyte* pDiag;

                    if (diag == MatrixDiagonalMode.NonUnitDiagonal)
                    {
                        sbyte tempuri = (sbyte) 'N';
                        pDiag = &tempuri;
                    }
                    else
                    {
                        sbyte tempuri = (sbyte) 'U';
                        pDiag = &tempuri;
                    }

                    sbyte* pSide;

                    if (side == MatrixOperationSide.Left)
                    {
                        sbyte tempuri = (sbyte) 'L';
                        pSide = &tempuri;
                    }
                    else
                    {
                        sbyte tempuri = (sbyte) 'R';
                        pSide = &tempuri;
                    }

                    ManagedBlasProvider.dtrmm_(pSide, mode, pTrans, pDiag, &m, &n, &alpha, pA, &lda, pB, &ldb);
                }
            }
        }

        /// <summary>
        /// Solution of a triangular linear system with multiple right-hand sides. 
        /// </summary>
        /// <param name="side">Specifies on which side the triangular matrix a is to be multiplied.</param>
        /// <param name="triangleMode">Specifies whether the elements of the matrix a are stored in the upper or lower triangular part.</param>
        /// <param name="transA">Specifies the operation to be performed on the matrix a.</param>
        /// <param name="diag">Specifies whether or not a is unit triangular.</param>
        /// <param name="m">The number of rows in the matrix a transformed as specified by transA, and the matrix c.</param>
        /// <param name="n">The number of columns in the matrix b transformed as specified by transB, and the matrix c.</param>
        /// <param name="alpha">The scalar used to multiply the matrix-vector product.</param>
        /// <param name="a">Reference to the first element in a one-dimensional array containing the elements of the first matrix.</param>
        /// <param name="lda">Leading dimension of a.</param>
        /// <param name="b">Reference to the first element in a one-dimensional array containing the elements of the second matrix.</param>
        /// <param name="ldb">Leading dimension of b.</param>
        public static unsafe void Dtrsm(MatrixOperationSide side, MatrixTriangleMode triangleMode,
                                        TransposeOperation transA,
                                        MatrixDiagonalMode diag, int m, int n, double alpha, double[] a, int lda,
                                        double[] b, int ldb)
        {
            fixed (double* pA = a)
            {
                fixed (double* pB = b)
                {
                    sbyte* pTrans;

                    if ((transA == TransposeOperation.None) || (transA == TransposeOperation.Conjugate))
                    {
                        sbyte tempuri = (sbyte) 'N';
                        pTrans = &tempuri;
                    }
                    else
                    {
                        sbyte tempuri = (sbyte) 'T';
                        pTrans = &tempuri;
                    }

                    sbyte* mode;

                    if (triangleMode == MatrixTriangleMode.Lower)
                    {
                        sbyte tempuri = (sbyte) 'L';
                        mode = &tempuri;
                    }
                    else
                    {
                        sbyte tempuri = (sbyte) 'U';
                        mode = &tempuri;
                    }

                    sbyte* pDiag;

                    if (diag == MatrixDiagonalMode.NonUnitDiagonal)
                    {
                        sbyte tempuri = (sbyte) 'N';
                        pDiag = &tempuri;
                    }
                    else
                    {
                        sbyte tempuri = (sbyte) 'U';
                        pDiag = &tempuri;
                    }

                    sbyte* pSide;

                    if (side == MatrixOperationSide.Left)
                    {
                        sbyte tempuri = (sbyte) 'L';
                        pSide = &tempuri;
                    }
                    else
                    {
                        sbyte tempuri = (sbyte) 'R';
                        pSide = &tempuri;
                    }

                    ManagedBlasProvider.dtrsm_(pSide, mode, pTrans, pDiag, &m, &n, &alpha, pA, &lda, pB, &ldb);
                }
            }
        }
    }
}